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1. INTRODUCTION 

Tliis report- gives a description of the work which lias been undertaken at. the University of Wales Swansea 
under NASA Research Grant NAGW — 2962. The objective of the research program was to produce an initial 
capability for the efficient simulation of 3D transient inviscicl compressible flows on unstructured tetrahedral 
meshes. There was also an objective of investigating the simulation of flows involving oscillating boundaries. 


2. WORK ACCOMPLISHED 

Over the past ten years, the views of the aerospace industry on the use of unstructured grid methods for the 
simulation of aerodynamic flows have changed considerably. In the early eighties, these methods were greeted 
with a large degree of suspicion. At that. time, the main causes for concern were the perceived difficulty of 
generating unstructured grids of high quality and the question of the accuracy of the results of computations 
made on general unstructured grids. In the intervening period, researchers have made major strides in the 
areas of grid generation, accuracy and adaptivity so that, nowadays, most major aerospace companies have, 
or are working towards the acquisition of. an unstructured grid aerodynamics capability. In addition, it is 
already recognised by industrial users that the unstructured approach will have an important, role to play 
in the development of any successful multi-disciplinarv (i.e. fluid dynamics, electro-magnetics, structural 
mechanics, optimisation) analysis capability (Shankar. Hall. Mohammadian and Rowell 1993). 

We will begin by oulining a general framework for the development, of algorithms for the solution of the 
inviscid flow equations on general unstructured tetrahedral grids. The approach which will be followed 
employs an edge based representation of the grid. It. will be shown that., as well as being a convenient 
mechanism for algorithm development., the edge based representation has significant, implications to the 
computational efficiency of the flow solvers, both in terms of computer time and memory requirements. 
Practical algorithms which have been utilised to date for the simulation of steady flows are generally based 
upon the use of explicit time stepping methods. In the structured grid context, more efficient methods 
have been developed which are based upon explicit time stepping with multigrid acceleration (Jameson 
1983). However, the implementation of the multigrid idea within the unstructured grid environment is not 
straightforward and certain difficulties need to be overcome. The use of explicit methods for the simulations 
of truly transient, flows can lead to significant computational requirements. The extension of the multigrid 
approach to handle transient problems is considered and the method is employed in the analysis of the flow 
over an oscillating wing. 

2.1 Equations of Inviscid Flow 

The equations of interest here are the 3D unsteady compressible Euler equations, expressed in a cartesian 
coordinate system Ox ixnx 3 in the conservation form 
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where the sumation convention is employed and where 
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is the unknown vector of conservation variables. The inviscid flux vectors are defined by 
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Here p, p and f denote t he pressure, density and specific total energy of the fluid respectively, while the 
component of the fluid velocity in the direction xj is denoted by uj . The equation set. is completed by the 
addition of the perfect gas equation of state 

/> = (-••- l)p(€-0.r)„ r ;,) (4) 


where 7 denotes the ratio of the specific heats. 

For the purpose of deriving algorithms for the solution of these equations in a region Q bounded by a closed 
surface T. using general unstructured tetrahedral grids, it is convenient to start, from the alternative weak 
variational formulation: find U(x,f) such that 

[ = [ f [ F (n, \vdr ( 5 ) 

Jn dt J n dxj ,/ r 

for every suitable weighting function \V(x) and every t > 0. subject to U(x.U) = Uo(x). In equation (5), 
the quantity F ( " 1 is defined In' 

F ( "’ = FV (6) 

where n = ( « j , no , /? 3 ) denotes the unit outward normal to T. The exact form of F (n| will depend on both 
the boundary condition being simulated and the local value of the solution. 

2.2 Unstructured Grid Solution Algorithms 

The domain Q is discretised using linear tetrahedral elements. In the standard finite element approach, a 
piecewise linear approximate solution U* is constructed in the form 

U*=U jNj(x) (7) 


where iVj(x) denotes the piecewise linear finite element shape function associated with node J of the grid 
and Uj = Uj(<) denotes the value of the approximate solution at node The nodal values are obtained 
from a standard Galerkin method, which is based upon the weak formulation of the problem in which 
the approximation is constructed by employing the shape functions themselves as the weighting functions 
(Johnson. 1987). Here, this is expressed as: find U* such that 

[ ^-N,dQ = f [ F (rt, .V/dr (8) 

Jn dt Jn dxj J r 

for every node I in the grid. The assumed form of U“ from equation (7) is inserted into this expression and 
the left hand side directly evaluated as 


L 


du m 

dt 


N,dQ = Mu 


dUj 

dt 


(9) 


where M/j denotes the entries in the standard consistent mass matrix. The right hand side terms in equation 
( 8 ) are evaluated by making use of the local nature of the finite element shape functions. The process is 
illustrated here for the domain integral only. Suppose that each edge in the grid is numbered and that the 
information about the grid is given in terms of the nodes which are joined by these edges. Suppose also that 
edge S in the grid joins the interior node / with node 5 t . If it is assumed that the fluxes are also linearly 
interpolated in terms of their nodal values i.e. 


F J ' « F J ( U j ) Nj ( x ) = F } j Nj ( x ) 
then it can lie readily shown (Peraire, Peiro Morgan 1993) that 
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where Cj ^ denotes the weight that must be applied to the sum of the edge nodal fluxes in order to obtain 
the contribution made to node [ by the edge which joins nodes / and 5j With an edge based representation 
for the grid, the discretised equation at. an interior node / can therefore he written as 


M ,j 


d\Jj 

dt 


( F '< + n.) 


( 12 ) 


It can be readily demonstrated from this expression that the flux at node / does not contribute to the 
discretised equation at. node I . so that the application of the Galerkin method leads to a scheme which is 
central difference in character. In addition, the total contribution made by any interior edge in the grid is 
zero, so that the scheme can immediately be seen to be conservative in this sense. 

To produce practical schemes for the simulation of compressible flows, we can employ the approach which is 
frequently adopted in ID and replace the actual fluxes in equation (12) by appropriate consistent numerical 
fluxes. To achieve this in a convenient manner, we begin by defining the vector 

C/.s, = (C}. Sl ,C}, Sl .C : ,\ Si ) ( 13 ) 

and the quantities 

£/.s, = I C/.5, I S}.s l =<?,.s t /Ci.s t (14) 

Then, the equation at node I may be expressed as 

= E £ ' < 15 ) 

Ssl 

where 

= F J ,Sj Si ^s t = HAs, ( 16 ) 

are the actual (luxes, in the direction of the vector C/s,, associated with the two nodes belonging to the 
edge S and 

$/.s, = (?[ + .Tys, ) ( 17 ) 

denotes the actual flux associated with the edge 5. Practical solution schemes for the compressible Euler 
equations on general unstructured meshes follow by replacing $ by a consistent <£. which is constructed on 
the basis of ID arguments. 

Centered Schemes 

Finite Element Lax-Wendroff Method 

This is a simple finite element scheme for the solution of the compressible Euler equations on unstructured 
grids and is the finite element equivalent, of the finite difference Lax-Wendroff method. For a problem 
involving only the spatial direction xq. the Lax-Wendroff flux function is defined as (Hirsch 1990) 


# = *'^ As(F ^‘" F}) 


( 18 ) 


where At is a time step, A* = dF l /d\J is the Jacobian matrix of F 1 and lig is the length of the edge 5. By 
direct extension to the multi-dimensional case, a Lax-Wendroff type scheme, for general tetrahedral grids, 
is written in the form 

MuAVj = A<££/. 5 ,*'; 5i ( 19 ) 

Ssl 


where 


In equation ( 19), 


* = $ - ^-a 5 (^ 5i -/■/) 

2 'IS 


AU = U" +l - U'‘ 


(20) 
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the superscript- n denotes the time level. A/ = /„ +t — and Ac - tlT /7/U. For the simulation of transient 
problems, one possible approach is to solve equation (19) by explicit iteration (Morgan Peraire 1987), 
though this will prove to be an expensive procedure in general. For steady flow simulations, the consistent 
mass matrix M is replaced by the diagonal lumped mass matrix Ml to produce a truly explicit scheme. 

.A Multi-Stage Method 

All approach which has proved to lie extremely successful for the simulation of compressible transonic Euler 
flows is the use of a central difference discretisation in space together with the explicit addition of a stabilising 
scalar high order added artificial viscosity (Jameson 1987). For a problem involving only the spatial direction 
X'i, a suitable flux function for a scheme of this form is 


&/.s, = */..s, - | -Vs 


U Sl - U s - I hs 


d U' 

(in 


s J 


(22) 


where the coordinate ij runs along the edge from node I to node Si- In equation ( 22 ), c /4 denotes a user- 
specified constant and | A | is the maximum eigenvalue, in absolute value, of the Jacobian matrix As- To 
implement, this procedure, nodal gradients are calculated, using equation ( 11 ), in the form 


[Ml],, 


du- 

dxj 


= -E c U( u ' + u *.) 

S:I 


(23) 


with j = 1. The required edge gradients are obtained as the average of the edge nodal gradients. Expression 
( 22 ) for the flux function can be shown to possess the property that the diffusion which is added is zero, on 
any grid, whenever U’ is linear. By direct, extension, we employ the flux function of equation (22) directly 
for the multi-dimensional case on general tetrahedral grids. 

The resulting equation can be advanced in time via the explicit multi-stage time stepping scheme 

U‘ 0) = U? 

U, fc 1 = uy + a t A t [Ml]7/ E ’ * = 1 (24) 

,Stl 

Un+l _ u<m) 


In an alterative implementation (Swanson & Turkel 1990). the scalar diffusion coefficient. | A 5 | , in equation 
( 22 ) is replaced by a matrix coefficient, leading to the numerical flux function 


‘fc/.Si = ‘fc/.Xi — d-1 


A 5 (U,.U 5l ) 


U Sl -U s - 



(25) 


where | As (U,,Us, ) | denotes the standard Roe matrix (Roe 1981) evaluated in the direction C/.s, • 
Shock Capturing 

Additional explicitly added low order artificial viscosity is needed if either of the above schemes is to be used 
successfully for the simulation of flows containing discontinuities. By comparison with the one dimensional 
case, it is apparent that the addition of diffusion in the multi-dimensional situation can be simply achieved 
by the use of the flux function 


4/.s, = i/,5, - d 2 *s I As I [Us, - U,] (26) 

where d 2 is a user defined constant. To maintain the accuracy of the original schemes in the smooth regions 
of the flow, an edge pressure switch 3 r 5 has been built- into the definition of this additional diffusion. This 
pressure switch is designed to be significant only near any flow discontinuities and to be very small in the 
smooth flow regions (Peraire. Morgan. Peiro & Zienkiewicz 1987). The alternative expression 

£/.s, = $/.s, - d?7r s | A 5 (U,,U 5 , ) | [U s , - U,] (27) 
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iii which tin? scalar diffusion coefficient, in equation (26) is replaced hy a matrix diffusion coefficient., may 
also be employed here. 

Higher Resolution via Flux Corrected Transport 

Flux corrected transport. (FCT) is a technique introduced originally (Boris k Book 1973. Zalesak 1979) 
to produce high resolution monotone solutions by combining the results of a low order monotone scheme 
with those of a higher order scheme. Unstructured grid implementations have proved successful (Parrott. k 
Christie 198(5. Morgan, Peraire k Lohner 1988). particularly for transient simulations. The approach only 
needs slight modification for implementation in the curent context. For example, for the Lax-YVendroff type 
scheme, let AU H be the basic (high order) solution increment and let AU 1 be the (low order) increment 
computed with the same algorithm but. with the addition of the unswitched (i.e. ;r<? - 1) shock capturing 
diffusion of equation (26). If this added diffusion is denoted by D. it follows that these increments are related 
by 

AU" = AU L - D (28) 

It is apparent, from equation (26) that the diffusion D is formed as the sum of edge contributions i.e. 

D = °/.5. (29) 

Sc I 

The idea behind FCT is to limit the amount of artificial diffusion which is removed from the low order 
solution to produce the high order solution in equation (28) by imposing the requirement that, the high order 
solution should remain monotone. This is achieved by firstly identifying maximum and minimum allowable 
values for the solution at. node /. We define 

U"=U"+AU W U L = U" + AU‘ (30) 

and the maximum and minimum allowable nodal values are determined as 

U' r nax = max (Vf , U£, , U? . U?, ) 

i£/ , . (31) 

Uy n = mn(Uf,US l ,U7.US l ) 

Edge limiters ti’s are determined so as to ensure that the increment computed according to 

AU = AU £ - Y, V’sD/.s, (32) 

Sc I 

is such that the new solution at. each node lies within the range specified by these maximum and minimum 
values. In practice, the edge limiter is determined from the requirement that, variables such as density, 
pressure and Mach number should remain monotone. As well as producing high resolution solutions, the 
FCT procedure acts as a stabilising influence on centered schemes when they are applied in the hypersonic 
range. Thus, with the addition of FCT, the explicit, multi-stage time stepping scheme has been employed 
to solve the problem of steady flow over a space shuttle type configuration at a Mach number of ten and at 
thirty degrees angle of attack. A view of the surface grid used in the simulat ion is shown in Figure la and the 
computed distribution of the Mach number contours can be seen in Figure lb. The complete unstructured 
grid, for one half of the configuration, consisted of 1,462,189 elements and 269.780 nodes. 

Upwind Methods 

By similar direct extension of schemes developed for the solution of the compressible Euler equations in ID, 
a first order upwind method for general unstructured grids can be constructed by the use of the flux function 

*/,s. = *r.s, - | AsfU/.UsJ | [U Sl -Ur] (33) 

Higher order methods can be derived by employing MUSCL (van Leer 1979) or TVD concepts (Harten 1983, 
1984). 
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ill the MUSC'L scheme (Vahdati. Morgan k Peraire 1992). the gradient, of the solution is computed at the 
nodes, as in equation (23). and the solution is linearly reconstructed for each edge. Slope limiting is applied 
to ensure that the reconstructed solution does not go above (below) the maximum (minimum) of the value 
of the solution at. the nodes of all elements associated with the two nodes of the edge. With a reconstructed 
solution, the scheme employs the numerical flux function 

*/.s, = (/■+ + /•-)- jA s (U+.U-)| [U^-U+] (34) 

where the superscripts + and - denote mid-edge states, which are determined from the reconstructed 
solution. 

In a TVD implementation (Lyra. Morgan. Peraire and Peiro 1993), points 5 ■, and S 3 are located in the grid 
such that (a) Si and 63 lie on the line of the edge S (b) St and S\ are equidistant, from / (c) Si lies on the 
opposite side of the point / to the point 5i (d) I and S 3 are equidistant from ,S'i (c) 5 3 lies on the opposite 
side of the point Si to the point I. The TVD flux is computed as (Yee 1989) 

§i,s t = $/.s, - nr 1 I A I Q (35) 

where the Roe matrix is decomposed according to 

| As | = R _1 | A | R (36) 

In this expression R denotes the matrix whose columns are the right eigenvectors of the Roe matrix and 
| A | is the diagonal matrix whose entries are the absolute values of the eigenvalues of A.s. I 11 equation (35), 

Q = Pi. Si - minmod (a’s,./, U’/.s, • <-'s,,s 3 ) (37) 

where 

Pi.Si = R(Us, - U/) (38) 

With either of these schemes, the solution can be advanced in time to steady state by any suitable time 
stepping procedure. 

Grid Data Structure 

For a general 3D grid of linear tetrahedral elements, the number of elements is approximately 5.5 times 
the number of nodes, the number of triangular faces is approximately twice the number of elements and 
the number of edges is approximately equal to the number of elements plus the number of nodes. Using 
these relations, it is possible to show that, in addition to being useful for the creation of flow algorithms, 
the adoption of the edge based representation for the unstructured grid has important computational 
implications. This is illustrated in the table below, which compares the storage and cpu requirements for an 
explicit 


Data Structure 

Storage/Node 

cpu /node 

Element 

25* 

1.0* 


71.5 

0.88 

Face 

58* 

1.4* 


88 

1.2 

Edge 

32.5 

0.6 


"With recomputation of the geometry 


centered unstructured grid solution algorithm for the solution of the compressible Euler equations. The 
grid data structures which are considered in this table are the traditional finite element data structure, the 
face data st ructure, which is frequently favoured by the developers of finite volume algorithms, and the edge 
data structure, which has been considered in detail above. It can be observed that the computer code based 
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upon the edge data, structure for the grid is the most efficient in terms ot cpu requirements and is close to 
being optimal in terms of its storage demands. 

2.3 Convergence Acceleration 

The solution of complex steady 3D flows on unstructured grids is most commonly achieved by explicit time 
stepping. Although this approach is robust, it can prove expensive in terms of cpu requirements, as a large 
number of time steps is frequently needed to produce a converged solution. There is therefore the requirement 
to search for more effective solution approaches, especially with the perceived increasing industrial demand 
for the development of a capability for the solution of large scale viscous (low problems. 

Multigrid Methods 

Multigrid acceleration coupled with an explicit multi-stage time stepping scheme has proved to be an 
extremely effective method for the solution of the compressible Euler equations on structured grids (Jameson 
1983 ) . The basic idea of the multigrid method is to use the residuals computed on a given grid to drive a time 
marching scheme on a coarse grid. The corrections to the unknowns computed on the coarser grid are then 
added to the fine grid solution. The time marching scheme on the coarse grid may itself be accelerated by 
the use of an even coarser grid and. in this wav, the multigrid concept can be extended to incorporated any 
number of coarse grids. Advancing the solution on the coarse grids is computationally relatively inexpensive, 
because of the reduced number of unknowns and the larger allowable time step size. The necessary int.ergrid 
transfers are readily accomplished because of the nested nature of the grid sequence. 

The development of appropriate intergrid transfer operators is a complication that has to be addressed when 
attempting to implement a similar multigrid acceleration scheme within the unstructured grid environment 
(Mavriplis 1988. Mavriplis 1991. Peraire, Morgan and Peiro 1993). 

2.4 Transient Problems 

Although the time dependent form of the equations has been employed, the main thrust of the arguments 
presented above has been directed at the simulation of steady flows, where the solution is obtained via a 
false transient. The extension of these techniques to the simulation of truly time dependent flows is now 
considered. The particular problem of simulating the flow over an oscillating wing has been chosen as an 
illustrative example. 

Flow Over an Oscillating Wing 

Two dimensional simulations of t ransonic inviscid flows over oscillating bodies have been performed using the 
full range of techniques currently available to the unstructured grid researcher. To allow for the movement of 
the boundary and the transient development of the solution, adaptive remeshing can be utilised. To improve 
the efficiency of the explicit computational procedure, on meshes exhibiting large variations in element size, 
a domain decomposition technique can be included. In the domain decomposition method, different time 
steps are employed in different areas of the grid, while maintaining the time accuracy of the overall solution 
procedure. Typical results are shown in Figures 2 and 3 for a standard AGARD test case involving flow at 
a freestream Mach number of 0.755 over an oscillating NACA0012 aerofoil. The procedure can also be used 
in a predictive mode, as shown in Figure 4 which displays the results of the flow over a NACA0012 aerofoil 
in which the rear portion only of the aerofoil is oscillating. 

Such comprehensive 2D computations place large demands on current computational resources. For 3D 
simulations, we therefore decided to investigate if we could produce an approach which is more in line with 
the current computational resources, at the expense of a certain degree of generality. To achieve this, we 
restricted consideration to simulations in which a moving boundary would actually be simulated by keeping 
the boundary and the mesh fixed and then applying a suitably modified boundary condition. This avoids 
the necessity of expensive remeshing procedures to allow for the boundary motion. To produce a further 
improvement in the computational performance, the explicit timestepping scheme is replaced by an implicit 
procedure. This allows for the adoption of larger timesteps, while still maintaining sufficient time accuracy. 
The solution of an implicit equation system resulting from the 3D Euler equations is a non-trivial task and 
the approach followed here has been to achieve this by the application of the multigrid method. 
Approximation to the Wing Boundary Condition 

To be able to impose the boundary condition at any point on a moving solid surface, we need to know the 
normal vector n to the surface and the normal velocity v„ of the surface at the point. The surface boundary 
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condition is then imposed by replacing the computed fluid velocity V by V according to 


V* = V — (V • n) n + (39) 

tangential normal 

The intention here is to keep the boundary fixed and to simulate its actual motion through a boundary 
condition of the form of equation (39). Suppose the angle of attack, o. of the wing is actually changing in 
time according to 

O' = oo sin(a;< + o) (40) 

The normal vector n at time t is computed from the original normal, tin. at time t = 0 as 

n = T(o)n 0 (41) 


where T(o) is a standard rotation matrix. The normal velocity of the wing surface is calculated from the 
total velocity Vr as 

v„ = W ■ n (42) 


with the total velocity VY determined as 


da 

V ’ T = r x Tt 


(43) 


where r denotes the position vector from the line of rotation to the point, on the surface in question. 

To investigate the usefulness of this form of the boundary condition, an initial investigation was undertaken 
in which a wing at an angle of attack was simulated (a) by inclining the flow direction and applying the 
standard inviscid boundary condition and (b) by a wing at zero angle of attack and applying the boundary 
condition of equations (39-43). The results show that the two approaches produce very similar results for 
reasonable angles of attack. For example, in Figure 5. we compare the pressure distributions, computed by 
the two approaches, over a section of an ONERA M6 wing for an angle of attack of four degrees. Results of 
this form give us confidence that the proposed approach can be used to give meaningful results for transient 
simulations in which the variation in the angle of attack is of this order. 

An Implicit Timestepping Scheme 

The spatial discretisation methods for the Euler equations which have been considered above lead to a 
discretised equation svstem 

dU 

— = R(U) (44) 

In the present context, we propose the use of an implicit three-level scheme in time (Jameson 1992) which 
can be expressed as 

3 U n+1 - 4 U" + U" -1 

^ =R(U"+‘) (45) 

to advance the solution. Equation (45) can be re-arranged and written in the form 


0 = R- (U" +l ) +S (U".U"-‘) 


(46) 


where 


and 


R- (U" + l ) = R(U" +1 ) - J3^U” +1 


S (UMJ’- 1 ) = 


2 U n U 


At 




2 At 


(47) 


(48) 


Equation (46) can now be solved by using an explicit multi-stage time stepping scheme to obtain the steady 
state solution of the equation 

dU ,,+ 1 

— = R- (U n + 1 ) +S(U'‘.U”- 1 ) 
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(49) 


with multigrid acceleration again employed to accelerate the convergence of the explicit time stepping. 

The numerical performance of the proposed procedure is demonstrated by considering the transient 
development of the flow over an oscillating ONERA M6 wing. In the computation, the wing is held fixed and 
its sinusoidal oscillation is simulated by modifying the boundary condition at the wing surface, according to 
equations (39-43). As we have indicated previously, this approach can be expected to be valid for oscillations 
of small amplit ude. Here the simulation involved an amplitude of four degrees and 40 multigrid cycles were 
used per timestep with 24 timesteps taken per cycle. A sequence of five meshes was employed in the multigrid 
procedure, with the form of the surface discretisations employed being indicated in Figure 6. The distribution 
of the pressure on the upper surface of the wing during the initial portion of the cycle is shown in Figure 7. 
The distribution of the pressure on the lower surface of the wing during the remainder of the cycle is shown 
in Figure 8. 


3. CONCLUSIONS 

A description of the research work undertaken under NASA Research Grant NAGW-2962 has been given. 
Basic algorithmic development work, undertaken for the simulation of steady 3D inviscid flows, has been 
used as the basis for the construction of a procedure for the simulation of truly transient flows in 3D. To 
produce a viable procedure for implementation on the current generation of computers, moving boundary 
components are simulated by fixed boundaries plus a suitably modified boundary condition. Computational 
efficiency is increased by the use of an implicit timestepping scheme in which the equation system is solved 
by explicit multistage timestepping with multigrid acceleration. The viablity of the proposed approach has 
been demonstrated by considering the application of the procedure to simulation of a transonic flow over an 
oscillating ONERA M6 wing. 
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Figure 3 

The simulation of inviscid flow over an oscillating NACA0012. Free stream Mach number = 0.755. 
Distribution of the computed pressure coefficient over the aerofoil surface showing a comparison with 
experiment and with the results of a simulation in which the mesh was adapted by moving only. 
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Figure 4 

The simulation of inviscid flow over a NACA0012, with the rear portion of the aerofoil oscillating sinusoidally.. 
Free stream Mach number = 0.8. Details of the adapted meshes and computed pressure coefficient contours 
at different stages of the cycle. 
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Figure 5 

The simulation of inviscid transonic flow over an ONERA M6 wing at 4 degrees angle of attack. Comparison 
of the computed pressure distributions over a section obtained using an inclined stream (de) and the modified 
boundary condition (sm). 
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Figure 6 

The simulation of inviscid transonic flow over a sinusoidally oscillating ONERA M6 wing. The surface 
discretisation on (a) The finest grid, (b) An intermediate grid, (c) The coarsest grid. 
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